########################
# 16_mexico_pressure.R #
########################

#######################################################
### Pressure model, 20 chains on 10 samples, Mexico ###
#######################################################


rm(list=ls())

library (runjags)
library (coda)
library (reshape2)
library (random)
library (scales)


print(version)

# platform       x86_64-pc-linux-gnu         
# arch           x86_64                      
# os             linux-gnu                   
# system         x86_64, linux-gnu           
# status                                     
# major          4                           
# minor          3.0                         
# year           2023                        
# month          04                          
# day            21                          
# svn rev        84292                       
# language       R                           
# version.string R version 4.3.0 (2023-04-21)
# nickname       Already Tomorrow            


# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)


all.seeds <- matrix (data= c(9198, 2430, 7696, 8307, 5173, 1881, 2233, 9302,
                             4049, 5652, 1938, 3437, 7129, 7921, 1140, 4605, 
                             8003, 1946, 9408, 5692),
                     ncol=2, byrow = FALSE)

row <- c(1:10)


for (j in 1:2){
  for (i in 1:10){
    which.sample <- row[i] # for sample_mean, use "mean"
    which.process <- ifelse (j==1, "a", "b")
    this.seed <- all.seeds[which.sample, j]
    # for sample_mean, use this.seed <- 1971 for a or 1972 for b
    
  load (paste0 ("mexico_rollcall_long_format_", which.sample, ".RData")) 

  load ("mexico_initial_values.Rda")

  potential.fix.bill <- c(22,89,144,150,174, 217)

  pressure<-  function(x1, x2, x3) { 
    k<- rep(0, length(x1))
    for (i in 1: length(x1)){
      k[i]<-ifelse(max(x1[i], x2[i], x3[i])>0, max(x1[i], x2[i], x3[i]), min(x1[i], x2[i], x3[i]))
    }
    k
  }

  pressure.data.vector <- dump.format(list(y=good.RC$RC
                                  , fix.bill=potential.fix.bill
                                  , n.legs=max(good.comp.party$final.leg.no)
                                  , n.item=max(as.numeric(good.comp.party$vote))
                                  , n.obs=nrow(good.RC)
                                  , obs=seq(1:nrow(good.RC))
                                  , vote=as.numeric(good.comp.party$vote)
                                  , dep=good.comp.party$final.leg.no
                , P1=good.R7$pressure #PRI
                , P2=good.R8$pressure #PAN
                , P3=good.R9$pressure #PRD
                , P4=pressure(good.R1$pressure, good.R2$pressure, good.R3$pressure) #Interaction with PR
                , P5=good.R67$pressure   #PR
                , P6=pressure(good.R27$pressure, good.R26$pressure, good.R28$pressure)
                , P7=good.R64$pressure  #poverty dummy
                , P8=pressure(good.R69$pressure, good.R68$pressure, good.R70$pressure)
                , P9=good.R65$pressure    # year
                , P10=pressure(good.R30$pressure,good.R29$pressure, good.R31$pressure)
                , P11=good.R75$pressure   #victory margin dummy
                , P12=pressure(good.R51$pressure, good.R50$pressure,good.R52$pressure)
                , P13=good.R72$pressure  #sta.muni.leg
                , P14=pressure(good.R48$pressure, good.R47$pressure, good.R49$pressure)
                , P15=good.R77$pressure  #CabinetExp
                , P16=pressure(good.R37$pressure, good.R33$pressure, good.R41$pressure)
                , P17=good.R71$pressure   #prezposition
                , P18=pressure(good.R45$pressure, good.R44$pressure, good.R46$pressure)
                , P19=good.R73$pressure    #commPres
                , P20=pressure(good.R36$pressure, good.R32$pressure,good.R40$pressure)
                , P21=good.R56$pressure    #pcbrequest
  ))


  pressure.parameters = c("theta", "delta", "alpha", "beta", "eta", "lambda", "deviance")

  pressure.only.v="model { 
     for (i in 1:n.obs) {
         y[i] ~ dbern(pi[i])
         probit(pi[i]) <- beta[vote[i]]*theta[dep[i]] +eta[vote[i]]*delta[dep[i]]- alpha[vote[i]] + lambda[1]*P1[obs[i]] +
                        lambda[2]*P2[obs[i]] + lambda[3]*P3[obs[i]] + lambda[4]*P4[obs[i]] +
                        lambda[5]*P5[obs[i]] + lambda[6]*P6[obs[i]] + lambda[7]*P7[obs[i]] +
                        lambda[8]*P8[obs[i]] + lambda[9]*P9[obs[i]] + lambda[10]*P10[obs[i]] + 
                        lambda[11]*P11[obs[i]] + lambda[12]*P12[obs[i]] + lambda[13]*P13[obs[i]] + 
                        lambda[14]*P14[obs[i]] + lambda[15]*P15[obs[i]] + lambda[16]*P16[obs[i]] + 
                        lambda[17]*P17[obs[i]] + lambda[18]*P18[obs[i]] + lambda[19]*P19[obs[i]] +
                        lambda[20]*P20[obs[i]] + lambda[21]*P21[obs[i]]
   } 
   # PRIORS
   for(j in 1:21) { lambda[j] ~ dnorm(0,1) } 
  # PRIORS
  for (j in 1:n.item){ alpha[j] ~ dnorm(0,1) }   

  # Alternative identification: Beta (discrimination, dimension 1)
  for(j in 1:(fix.bill[1]-1)) { beta[j]  ~ dnorm(0,1) }
  beta[fix.bill[1]]~ dnorm(0,1)
  for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { beta[j]  ~ dnorm(0,1) }
  beta[fix.bill[2]] <- 1.5
  for(j in (fix.bill[2]+1):(fix.bill[3]-1)) { beta[j]  ~ dnorm(0,1) }
  beta[fix.bill[3]] ~ dnorm(0,1)  #<- 2
  for(j in (fix.bill[3]+1):(fix.bill[4]-1)) { beta[j]  ~ dnorm(0,1) }
  beta[fix.bill[4]] <- -1
  for(j in (fix.bill[4]+1):(fix.bill[5]-1)) { beta[j]  ~ dnorm(0,1) }
  beta[fix.bill[5]] <-  2.5
  for(j in (fix.bill[5]+1):(fix.bill[6]-1)) { beta[j]  ~ dnorm(0,1) }
  beta[fix.bill[6]] ~ dnorm(0,1)
  for(j in (fix.bill[6]+1):n.item) { beta[j]  ~ dnorm(0,1) }
  # Eta (discrimination, dimension 2)
  for(j in 1:(fix.bill[1]-1)) { eta[j]  ~ dnorm(0,1) }
  eta[fix.bill[1]] ~ dnorm(0,1)
  for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { eta[j]  ~ dnorm(0,1) }
  eta[fix.bill[2]] <- 2
  for(j in (fix.bill[2]+1):(fix.bill[3]-1)) { eta[j]  ~ dnorm(0,1) }
  eta[fix.bill[3]] ~ dnorm(0,1) #
  for(j in (fix.bill[3]+1):(fix.bill[4]-1)) { eta[j]  ~ dnorm(0,1) }
  eta[fix.bill[4]] <- -2 #
  for(j in (fix.bill[4]+1):(fix.bill[5]-1)) { eta[j]  ~ dnorm(0,1) }
  eta[fix.bill[5]] <- 0 #
  for(j in (fix.bill[5]+1):(fix.bill[6]-1)) { eta[j]  ~ dnorm(0,1) }
  eta[fix.bill[6]] ~ dnorm(0,1) #
  for(j in (fix.bill[6]+1):n.item) { eta[j]  ~ dnorm(0,1) }

  # ideal points
  for(i in 1:n.legs)  { theta[i] ~ dnorm(0,1) }
  for(i in 1:n.legs)  { delta[i] ~ dnorm(0,1) }
  ## Notes
  # In this model we only include pressure, but make no attempt to model abstentions
  }"


  ##Initials for identification
  set.seed (this.seed)
  eta1 = initials$eta + runif(length(initials$eta), -0.2, 0.2)
  beta1 = initials$beta + runif(length(initials$beta), -0.2, 0.2)
  eta1[potential.fix.bill[4:5]] <- c(NA,NA)
  eta1[potential.fix.bill[2]] <- c(NA)
  beta1[potential.fix.bill[2]] <- c(NA)
  beta1[potential.fix.bill[4:5]] <- c(NA,NA)

  pressure.inits.1 <- function (){
    dump.format (
      list(
        theta = initials$theta + runif(length(initials$theta), -0.2, 0.2),
        delta = initials$delta + runif(length(initials$delta), -0.2, 0.2),
        alpha = initials$alpha + runif(length(initials$alpha), -0.2, 0.2),
        eta = eta1,
        beta = beta1,
        lambda = rnorm(21,0,1)
        ,'.RNG.name'="base::Wichmann-Hill"
        ,'.RNG.seed'= this.seed
      )
    )
  }

  date()
  set.seed (this.seed)

  pressure.model <- run.jags(
   model=pressure.only.v,  
   monitor=pressure.parameters,
   data=pressure.data.vector,
   inits=list(pressure.inits.1()),
   thin=30, burnin=4000, sample=100, 
   check.conv=FALSE, plots=FALSE)

   chainsPress <- mcmc.list(list (pressure.model$mcmc[[1]]))


  date()

  save.name <-  paste0 ("mexico_pressure_final_result_sample_", which.sample,
                        which.process, ".Rda")

  chainsPress <- pressure.model$mcmc[[1]]  
  # Uncomment next line to save files
  # save (chainsPress, file=save.name)
}
}




